
## load relevant libraries 
library(devtools)
install_github("naoki-egami/factorEx", dependencies=TRUE)
#install.packages("factorEx")
library(factorEx)


#load working data from main analysis

setwd("~/Dropbox/women_compaign_project/Paper1/jop_submission/revised_paper/replication_files/")
master_ps <- read.csv("~/Dropbox/women_compaign_project/Paper1/jop_submission/revised_paper/replication_files/master_ps.csv")

names(master_ps)


# marginals 


gender_dist_marginals <- list(gender=c("Male"=0.76, "Female"=0.24),party=c("Independent"=0.56, "Minor"= 0.26, "Major"=0.26), promises= c("Boreholes"=0.33, "Education"=0.33, "Roads"= 0.33), education=c("Secondary"=0.55, "University"=0.45), profession=c("Teacher"=0.42, "Maize farmer"=0.05, "Major business owner"=0.53))



## using the factEx package to include marginals

sudat <- master_ps[,c("outcome_binary_resp","gender","party","promises","education","profession","KEY","PARENT_KEY","tr_video_lv")]

sudat <-na.omit(sudat)


# make factors 

sudat$gender <- factor(sudat$gender,levels = c("Male", "Female"))
sudat$party <- as.factor(sudat$party)
sudat$promises <- as.factor(sudat$promises)
sudat$education <- as.factor(sudat$education)
sudat$profession <- as.factor(sudat$profession)
sudat$pair_id <- as.integer(substr(sudat$KEY, nchar(sudat$KEY) - 1, nchar(sudat$KEY) - 1))


sudat <- sudat[order( sudat[,8], sudat[,10] ), ]

sudat$pair_id_rv <- as.integer(rep(1:13434, each=2))



## Estimate gender effect in the full sample
set.seed(12970)

out_model <- 
  factorEx::model_pAMCE(formula = outcome_binary_resp ~ gender + party + promises + education + profession,
                        data = sudat, 
                        reg =  TRUE,
                        pair_id = sudat$pair_id_rv,
                        cluster_id = sudat$PARENT_KEY,
                        target_dist  = gender_dist_marginals, 
                        target_type = "marginal",
                        boot = 2000)



#summary(out_model, factor_name = c("gender"), sample = TRUE)

#plot(out_model, factor_name = c("gender"), diagnose = TRUE)


#estimating the pAMCE by treatment 

table(sudat$tr_video_lv)

subdat_control <- sudat[sudat$tr_video_lv=="Control",]

subdat_progress <- sudat[sudat$tr_video_lv=="High viability",]

subdat_discr <- sudat[sudat$tr_video_lv=="Low viability",]


out_model_control <- 
  factorEx::model_pAMCE(formula = outcome_binary_resp ~ gender + party + promises + education + profession,
                        data = subdat_control, 
                        reg =  TRUE,
                        pair_id = subdat_control$pair_id_rv,
                        cluster_id = subdat_control$PARENT_KEY,
                        target_dist  = gender_dist_marginals, 
                        target_type = "marginal",
                        boot = 2000)

#summary(out_model_control, factor_name = c("gender"), sample = TRUE)


#plot(out_model_control, factor_name = c("gender"), diagnose = TRUE)

out_model_progress <- 
  factorEx::model_pAMCE(formula = outcome_binary_resp ~ gender + party + promises + education + profession,
                        data = subdat_progress, 
                        reg =  TRUE,
                        pair_id = subdat_progress$pair_id_rv,
                        cluster_id = subdat_progress$PARENT_KEY,
                        target_dist  = gender_dist_marginals, 
                        target_type = "marginal",
                        boot = 2000)

#summary(out_model_progress, factor_name = c("gender"), sample = TRUE)

out_model_d <- 
  factorEx::model_pAMCE(formula = outcome_binary_resp ~ gender + party + promises + education + profession,
                        data = subdat_discr, 
                        reg =  TRUE,
                        pair_id = subdat_discr$pair_id_rv,
                        cluster_id = subdat_discr$PARENT_KEY,
                        target_dist  = gender_dist_marginals, 
                        target_type = "marginal",
                        boot = 2000)

#summary(out_model_d, factor_name = c("gender"), sample = TRUE)



## put results into table
restab <- rbind(summary(out_model, factor_name = c("gender"), sample = TRUE), summary(out_model_control, factor_name = c("gender"), sample = TRUE), summary(out_model_progress, factor_name = c("gender"), sample = TRUE),summary(out_model_d, factor_name = c("gender"), sample = TRUE))


sink("tab_pamce.tex")

stargazer::stargazer(restab,summary = FALSE)

sink()